close all 
clear all


%---------------------- Plech Data -------------------------

%s52=load('52nm.txt');
%s94=load('94nm.txt');
%f52=s52(:,1);da52=s52(:,2);
%f94=s94(:,1);da94=s94(:,2);
%T52=zeros(length(f52),1);
%T94=zeros(length(f94),1);

%tp=100e-15*(2/log(2))^(1/2);
%G=105e6;
%g=2e16;
%To=300;
%Ceo=70*To;

%Q52=1.731;
%Q94=1.956;
%r52=52/2*10^(-9);
%r94=94/2*10^(-9);

%g52=g*tp/Ceo;
%g94=g*tp/Ceo;
%b52=Q52*pi*r52^2/(Ceo*4/3*pi*r52^3*To).*f52;
%b94=Q94*pi*r94^2/(Ceo*4/3*pi*r94^3*To).*f94;
%e52=3*G/(r52*g);
%e94=3*G/(r94*g);

%Tguess=To;
%dT=0.25;
%for ii=1:length(T52)
%    check=10;
%    T=Tguess;
%    while check>=0.5
%        T=T+dT;
%        a52=(12.00269+0.00953*T-8.4*10^(-6)*T^2+5.43*10^(-9)*T^3)*10^(-6);
%        T52(ii)=da52(ii)./a52+To;
%        check=(T52(ii)-T)/T52(ii)*100;
%    end
%end

%for ii=1:length(T94)
%    check=10;
%    T=Tguess;
%    while check>=0.5
%        T=T+dT;
%        a94=(12.00269+0.00953*T-8.4*10^(-6)*T^2+5.43*10^(-9)*T^3)*10^(-6);
%        T94(ii)=da94(ii)./a94+To;
%        check=(T94(ii)-T)/T94(ii)*100;
%    end
%end

%theta52=(T52-To)./((b52./g52+1).*To-To);
%theta94=(T94-To)./((b94./g94+1).*To-To);

%atheta52=mean(theta52);
%atheta94=mean(theta94);

%stheta52=std(theta52);
%stheta94=std(theta94);

%---------------------------------------------------------------------------

vgamma_temp=load('sphere_10_15//vgamma.txt');
vbeta_temp=load('sphere_10_15//vbeta.txt');
veps_temp=load('sphere_10_15//veps.txt');

te_temp=load('sphere_10_15//e_eps.txt');
tp_temp=load('sphere_10_15//p_eps.txt');
tw_temp=load('sphere_10_15//w_eps.txt');

[yy,xx]=size(te_temp);
yye=yy;
te=zeros(yy*xx,1);
ee=te;eg=te;eb=te;
for ii=1:xx
    for jj=1:yy
        loc=(ii-1)*yy+jj;
        te(loc,1)=te_temp(jj,ii);
        ee(loc,1)=veps_temp(ii);
        eb(loc,1)=vbeta_temp(jj);
    end
end

[yy,xx]=size(tp_temp);
tp=zeros(yy*xx,1);
pe=tp;pg=tp;pb=tp;
for ii=1:xx
    for jj=1:yy
        loc=(ii-1)*yy+jj;
        tp(loc,1)=tp_temp(jj,ii);
        pe(loc,1)=veps_temp(ii);
        pg(loc,1)=vgamma_temp(jj);
    end
end

[yy,xx]=size(tw_temp);
tw=zeros(yy*xx,1);
we=tw;wg=tw;wb=tw;
for ii=1:xx
    for jj=1:yy
        loc=(ii-1)*yy+jj;
        tw(loc,1)=tw_temp(jj,ii);
        we(loc,1)=veps_temp(ii);
        wg(loc,1)=vgamma_temp(jj);
    end
end

%CE=[-2.6982;1.0062;-0.1580;0.3171;2.5752;-0.2360];
%CP=[-7.1;1;7.1;1;1/3;1;7;0.5;-1;0.1;1;1;1];
%CW=[-8.5;-0.58;0.1;1;1/3;1;7;0.5;-1;0.1;1;1;1;1;-0.58];

CE=[-2.6982; 1.0063;-0.1580;0.3171;2.5753;-0.2360];
CP=[-7.3354;1.0366;6.3289;1.0000;0.2221;0.5810;10.9347;0.3590;-2.4720;
    3.5307;0.7242;-0.9397;1.9275];
CW=[-3.8988;0.1338;-0.1235;0.4859;0.6325;0.4670;8.1412;-0.0942;-3.5386;
    4.0336;-0.0002;1.0000;5.2826;5.2538;-0.8207];   

logte=log(te);
logtp=log(tp);
logtw=log(tw);
epara(:,3)=ee;
epara(:,2)=log(eg);
epara(:,1)=log(eb);
ppara(:,3)=pe;
ppara(:,2)=log(pg);
ppara(:,1)=log(pb);
wpara(:,3)=we;
wpara(:,2)=log(wg);
wpara(:,1)=log(wb);

options.MaxFunEvals=1*10^300;
options.TolX = 1e-30;
options.TolFun = 1e-30;

[c,res]=lsqcurvefit(@rTfunce2,CE,epara,logte);
cfite=c;
[c,res]=lsqcurvefit(@rTfuncp2,CP,ppara,logtp);
cfitp=c;
[c,res]=lsqcurvefit(@rTfuncw2,CW,wpara,logtw);
cfitw=c;

logte_out=rTfunce2(cfite,epara);
logtp_out=rTfuncp2(cfitp,ppara);
logtp_out2=rTfuncp2(CP,ppara);
logtw_out=rTfuncw2(cfitw,wpara);
logtw_out2=rTfuncw2(CW,wpara);

figure(1)
loglog(eb(1:yye),te(1:yye),'o');
hold on
loglog(eb(1:yye),exp(logte_out(1:yye)));
set(gca, 'xscale', 'log','yscale', 'log')
xlabel('\beta','FontSize',16);
ylabel('\theta_e','FontSize',16);

figure(2)
%loglog(pg,tp(:,1),'.');
loglog(pg(1:yy),tp(1:yy),'ro');
hold on
loglog(pg(yy+1:2*yy),tp(yy+1:2*yy),'go');
loglog(pg(2*yy+1:3*yy),tp(2*yy+1:3*yy),'bo');
loglog(pg(3*yy+1:4*yy),tp(3*yy+1:4*yy),'mo');
loglog(pg(4*yy+1:5*yy),tp(4*yy+1:5*yy),'co');

loglog(pg(5*yy+1:6*yy),tp(5*yy+1:6*yy),'rv');
loglog(pg(6*yy+1:7*yy),tp(6*yy+1:7*yy),'gv');
loglog(pg(7*yy+1:8*yy),tp(7*yy+1:8*yy),'bv');
loglog(pg(8*yy+1:9*yy),tp(8*yy+1:9*yy),'mv');
axis([10^(-3) 10^7 10^(-6) 10^3]);

loglog(pg(1:yy),exp(logtp_out(1:yy)),'r');
loglog(pg(yy+1:2*yy),exp(logtp_out(yy+1:2*yy)),'g');
loglog(pg(2*yy+1:3*yy),exp(logtp_out(2*yy+1:3*yy)),'b');
loglog(pg(3*yy+1:4*yy),exp(logtp_out(3*yy+1:4*yy)),'m');
loglog(pg(4*yy+1:5*yy),exp(logtp_out(4*yy+1:5*yy)),'c');
%errorbar(g52,atheta52,2*stheta52,'ro');
%errorbar(g94,atheta94,2*stheta94,'ko');
set(gca, 'xscale', 'log','yscale', 'log')
xlabel('\gamma','FontSize',16);
ylabel('\theta_p','FontSize',16);

figure(3)
loglog(wg(1:yy),tw(1:yy),'ro');
hold on
loglog(wg(yy+1:2*yy),tw(yy+1:2*yy),'go');
loglog(wg(2*yy+1:3*yy),tw(2*yy+1:3*yy),'bo');
loglog(wg(3*yy+1:4*yy),tw(3*yy+1:4*yy),'mo');
loglog(wg(4*yy+1:5*yy),tw(4*yy+1:5*yy),'co');

loglog(wg(5*yy+1:6*yy),tw(5*yy+1:6*yy),'rv');
loglog(wg(6*yy+1:7*yy),tw(6*yy+1:7*yy),'gv');
loglog(wg(7*yy+1:8*yy),tw(7*yy+1:8*yy),'bv');
loglog(wg(8*yy+1:9*yy),tw(8*yy+1:9*yy),'mv');
axis([10^(-3) 10^7 10^(-7) 10^3]);

loglog(wg(1:yy),exp(logtw_out(1:yy)),'r');
loglog(wg(yy+1:2*yy),exp(logtw_out(yy+1:2*yy)),'g');
loglog(wg(2*yy+1:3*yy),exp(logtw_out(2*yy+1:3*yy)),'b');
loglog(wg(3*yy+1:4*yy),exp(logtw_out(3*yy+1:4*yy)),'m');
loglog(wg(4*yy+1:5*yy),exp(logtw_out(4*yy+1:5*yy)),'c');
set(gca, 'xscale', 'log','yscale', 'log')
xlabel('\gamma','FontSize',16);
ylabel('\theta_w','FontSize',16);

%%%% 95% confidence interval standard deviations %%%%

twosw=2*(sum(sum(((exp(logtw)-exp(logtw_out))/exp(logtw)).^2))/length(logtw))^0.5*100
twosp=2*(sum(sum(((exp(logtp)-exp(logtp_out))/exp(logtp)).^2))/length(logtp))^0.5*100
twose=2*(sum(sum(((exp(logte)-exp(logte_out))/exp(logte)).^2))/length(logte))^0.5*100

